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Abstract. We discuss a new Monte Carlo algorithm for the simulation of complex 
fluids. This algorithm employs geometric operations to identify clusters of particles 
that can be moved in a rejection-free way. It is demonstrated that this geometric 
cluster algorithm (GCA) constitutes the continuum generalization of the Swendsen- 
Wang and Wolff cluster algorithms for spin systems. Because of its nonlocal nature, 
it is particularly well suited for the simulation of fluid systems containing particles 
of widely varying sizes. The efficiency improvement with respect to conventional 
simulation algorithms is a rapidly growing function of the size asymmetry between 
the constituents of the system. We study the cluster-size distribution for a Lennard- 
Jones fluid as a function of density and temperature and provide a comparison 
between the generalized GCA and the hard-core GCA for a size-asymmetric mixture 
with Yukawa- type couplings. 

1 Introduction and Motivation 

The presence of multiple time and length scales constitutes one of the major 
problems in computer simulations of matter. In simulations that faithfully 
capture the dynamic evolution of a system, the fastest particles in the system 
dictate the required time resolution. If different types of particles with widely 
varying diffusion rates are present, then the slower particles may be unable 
to explore the entire configuration space during the course of the simulation, 
leading to ergodicity problems. 

In the computational study of complex fluids, such as colloidal suspen- 
sions, this problem frequently occurs, since such systems typically contain 
particles of different sizes. A concrete example is the 'nanoparticle haloing' 
phenomenon discovered by Lewis and coworkers [1] . This experimental work 
deals with a new approach to the stabilization of suspensions of micron- 
sized spherical particles, which tend to aggregate under the influence of their 
mutual van der Waals attraction. Conventional approaches to prevent this 
gelation, such as charge stabilization (variation of the pH to alter the surface 
charge of the particles), lead to complications in certain applications, e.g., 
in the formation of colloidal crystals, where the electrostatic repulsion pre- 
vents the close packing of the particles. It was found that these complications 
can be avoided by generating an effective colloidal repulsion through the ad- 
dition of highly-charged nanometer-sized particles to the suspension of the 
(near-neutral) colloids. A concrete explanation for the underlying mechanism 
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leading to the effective repulsions is currently lacking. Evidently, a computa- 
tional approach to this problem must involve both the microspheres and the 
nanoparticles, which typically differ by a factor 100 in diameter. Tradition- 
ally, the effective pair interaction between colloids (potential of mean force) 
is then calculated by "integrating out" the smaller species, e.g., by simulat- 
ing a system containing two spheres at fixed separation embedded in a sea 
of smaller particles [2,3]. Here, we discuss a new simulation algorithm [4] 
that is capable of explicitly incorporating both species and computing equi- 
librium properties of the suspension, without suffering from the disparity 
in time scales. This algorithm can be viewed as a continuum version of the 
widely-used cluster algorithms for lattice spin models [5,6]. 

2 Cluster Monte Carlo Algorithms 

Equilibrium Monte Carlo methods are aimed at obtaining static thermody- 
namic and structural properties by generating system configurations accord- 
ing to the Boltzmann distribution. Nonphysical moves may be employed in 
the underlying Markov process, which - in principle - offers an exquisite 
way to overcome the presence of multiple time scales. Thus, Monte Carlo 
algorithms are the method of choice for, e.g., simulations of critical phenom- 
ena and phase transitions. In the context of size-asymmetric fluids, collective 
moves have been devised in which groups of particles are moved simulta- 
neously. Unless these groups are identified in a careful way, which typically 
involves knowledge about the physical properties of the system, the Monte 
Carlo step will entail a large energy change and consequently have only a very 
small acceptance rate. Thus, while this approach has been used quite suc- 
cessfully in specific situations [7-10], it typically involves a tunable parameter 
and cannot be generalized in a straightforward fashion. 

The situation is rather different for spin models, for which Swendsen and 
Wang (SW) [5] introduced a cluster algorithm in which groups of paral- 
lel spins are identified in a probabilistic manner, based upon the Fortuin- 
Kasteleyn mapping of the Potts model onto the random-cluster model [11]. 
These groups can subsequently be flipped independently. This rejection-free 
algorithm (every completed cluster can be flipped without an additional eval- 
uation of the resulting energy difference) suppresses critical slowing down. 
Accordingly, its generalization to off-lattice fluids has been a widely-pursued 
goal. However, the SW algorithm - as well as the even more efficient single- 
cluster variant introduced by Wolff [6] - relies on invariance of the Hamilto- 
nian under a spin-inversion operation. For a fluid, this translates into particle- 
hole symmetry, which is only obeyed for lattice gases. Consequently, efficient 
off-lattice cluster algorithms have only been designed for a small number 
of specific fluid models [12,13], and it is not clear how these methods can 
be generalized. Hard-sphere fluids constitute another special case. Here ev- 
ery configuration without particle overlaps has the same energy. Dress and 
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Krauth [14] devised a strategy to create such configurations by means of 
geometric transformations. Their approach is particularly advantageous in 
size-asymmetric mixtures, but cannot be applied to systems with other in- 
teractions without supplementing it with a costly acceptance criterion. Nev- 
ertheless, a remarkable feature of this geometric cluster algorithm is that 
it relies on the invariance of the Hamiltonian under geometric operations. 
Here, we exploit this property to formulate a cluster Monte Carlo algorithm 
that is applicable to arbitrary pair potentials without the imposition of an 
acceptance criterion. 

3 Generalized Geometric Cluster Algorithm 
3.1 Single-cluster variant 

The geometric cluster algorithm as formulated by Dress and Krauth starts 
from a configuration of particles, with periodic boundary conditions. This 
configuration is rotated around an arbitrarily chosen pivot. Groups of par- 
ticles that overlap between the original and the rotated configuration are 
exchanged between these configurations. In practice, it is more convenient 
to construct only a single cluster and to carry out a point reflection on the 
fly [15]. For each cluster a new pivot is chosen. This constitutes the counter- 
part of the Wolff algorithm for spin models [6] . In the presence of a general, 
isotropic pair potential V(r) the cluster construction proceeds as follows (cf. 
Fig. 1): 

1. Choose a random pivot point. 

2. Choose the first particle i at random and move it from its original posi- 
tion Yi to its new position rj, via a point reflection with respect to the 
pivot. 
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Fig. 1. Two-dimensional illustration of the interacting geometric cluster algorithm. 
Open and shaded circles label the particles before and after the geometrical oper- 
ation, respectively. The small filled disk denotes the pivot, a) Initial configuration; 
b) construction of a new cluster via point reflection of particles 1-3 with respect to 
the pivot; c) final configuration. 
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3. Identify all particles that interact with i in its original position or in 
its new position. These particles j are considered for a point reflection 
with respect to the pivot. This reflection is carried out with a probability 
Pij = max[l - exp(-/?Ajj),0], where Ay = V(|r- - rj|) - Vflr* - r^l) 
and (3 = 1/fceT. Note that py solely depends on the interaction strength 
between i and j. 

4. Repeat step 3 in an iterative fashion for each particle j that is added to 
the cluster. If j is moved with respect to the pivot, then all its interacting 
neighbors that have not yet been added to the cluster are considered for 
inclusion as well. The cluster construction is completed once all interact- 
ing neighbors have been considered. 

It is instructive to compare this prescription to the Wolff cluster algo- 
rithm, in which a cluster of parallel spins is grown from a randomly chosen 
initial spin. Parallel spins with a ferromagnetic coupling constant K are added 
to the cluster with a probability 1 — exp(— 2K). The energy difference between 
the parallel and the antiparallcl configuration indeed equals 2K. An antipar- 
allel spin is never added to the cluster, in accordance with the fact that such 
a pair will be in a state of lower energy if only the first spin is flipped. 

This algorithm is ergodic, since each particle can be moved over an ar- 
bitrarily small distance. Namely, there is a non-vanishing probability that 
a cluster consists of only a single particle and the pivot can be located ar- 
bitrarily close to the center of this particle. Detailed balance is proven by 
considering a configuration X which is transformed into a new configura- 
tion Y by means of a cluster move. The energy change Ey — Ex results 
from every particle that is not included in the cluster but interacts with one 
or more particles that are part of the cluster. Each "broken bond" has a 
probability 1—pk, so that the total probability of forming a cluster is given 



The prefactor C accounts for the probability of choosing a specific pivot 
and for the probability of creating a specific arrangement of bonds inside the 
cluster. The total set {k} of broken bonds can be divided into two subsets. The 
broken bonds I that lead to an increase A; in pair energy have a probability 
exp(— pAi), whereas the broken bonds m that lead to a decrease in pair energy 
have a probability equal to unity. Accordingly, the transition probability can 
be written as 



by 



T(X^Y) = C\{{l-p k ) . 



(1) 
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The reverse move, in which the configuration Y is transformed into the con- 
figuration X by moving a cluster that is constructed in the same way, has a 
probability 



T(Y^X) = C]l(l-p k ), 



(3) 



where the sum runs over the same set {k} of broken bonds as in Eq. (1), 
but the sign of all energy differences has been reversed (indicated by pk)- 
Accordingly, the sum over A; in Eq. (2) is replaced by the negative sum over 
the complementary set {m}, 



T(Y ^X) = Cexp 



(4) 



The point reflection is a self-inverse operation, so that detailed balance is 
obeyed without the need to impose an acceptance criterion: 



T(X -> Y) 
T(Y X) 



— exp 



k 



exp(-pE x ) 



(5) 



Since energy differences are taken into account on a pair-wise basis during 
the construction of the cluster, rather than as a total energy difference after 
the cluster has been completed, the cluster is constructed in such a way that 
large energy differences are avoided. The clusters are representative of the 
actual structure of the system. 



3.2 Multiple-cluster variant 

In order to demonstrate that the generalized geometric cluster algorithm 
(GCA) indeed constitutes the off-lattice counterpart of the SW and Wolff 
cluster algorithms, it is instructive to formulate a multiple-cluster variant. 
This formulation yields a full decomposition of an off-lattice fluid configura- 
tion into stochastically independent clusters. In the following, we demonstrate 
how it can be phrased as a natural extension of the single-cluster version. 

First, a cluster is constructed according to the Wolff version of the GCA, 
with the exception that the cluster is only identified; particles belonging to 
the cluster are marked but not actually moved. The chosen pivot will also 
be used for the construction of all subsequent clusters in this decomposition. 
These subsequent clusters are built just like the first cluster, except that 
particles that are already part of an earlier cluster will never be considered 
for a new cluster. Once each particle is part of a cluster the decomposition is 
completed and each cluster is moved with a probability /. 

Although all clusters except the first are built in a restricted fashion, 
every individual cluster is constructed according to the rules of the Wolff 
formulation of the GCA. The exclusion of particles that are already part of 
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another cluster simply reflects the fact that every bond should be considered 
only once. If a bond is broken during the construction of an earlier cluster it 
must not be re-established during the construction of a subsequent cluster. 

In order to establish that this prescription is a true equivalent of the SW 
algorithm, we prove that each cluster can be moved (reflected) independently 
while preserving detailed balance. If only a single cluster is actually moved, 
this essentially corresponds to the Wolff version of the GCA, since each clus- 
ter is built according to the GCA prescription. The same holds true if several 
clusters are moved and no interactions are present between particles that 
belong to different clusters (the hard-sphere algorithm is a particular realiza- 
tion of this situation). If two or more clusters are moved and broken bonds 
exist between these clusters, i.e., a non- vanishing interaction exists between 
particles that belong to disparate (moving) clusters, then the shared broken 
bonds are actually preserved and the proof of detailed balance provided in 
the previous section no longer applies in its original form. However, since 
these bonds are identical in the forward and reverse move, the corresponding 
factors cancel out. This is illustrated for the situation of two clusters whose 
construction involves, respectively, two sets of broken bonds {fci} and {fe}. 
Each set comprises bonds I that lead to an increase in pair energy and bonds 
in that lead to a decrease in pair energy. We further subdivide these sets into 
external bonds that connect cluster 1 or 2 with the remainder of the system 
and joint bonds that connect cluster 1 and 2. Accordingly, the probability of 
creating cluster 1 is given by 

ci H(i-p i )=c 1 n a-Pi) n a-^)- ( 6 ) 

*e{'i} te{q^} je{l[ oint } 

Upon construction of the first cluster, the creation of the second cluster has 
a probability 

C2 n (i-w). (7) 

since all joint bonds in {Zj° mt } already have been broken. The factors C\ 
and Ci account for the probability of realizing a particular arrangement of 
internal bonds in clusters 1 and 2, respectively. Hence, the total transition 
probability for moving both clusters (upon fixing the pivot) is given by 



T\i = C\Ci exp 



ie{if«} jeljr} ne{;i oint } 



(8) 



In the reverse move, the energy differences for all external broken bonds have 
changed sign, but the energy differences for the joint bonds connecting cluster 
1 and 2 are the same as in the forward move. Thus, cluster 1 is created with 
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probability 



ci n c 1 -**) n (i-w) 

ie{mft} je{if nt } 
ie{mft} je{if nt } 



(9) 



where the p reflects the sign change compared to the forward move and the 
product over the external bonds involves the complement of the set {l° xt }. 
The creation probability for the second cluster is 



C2 n ( i -po=c 2 n cx p[+^] 

ie{m| xt } ie{m| xt } 

and the total transition probability for the reverse move is 



(10) 



T12 = CiC 2 exp 



+(3 ]T 2 



ie{m° xt } ie{m| xt } 

Accordingly, detailed balance is still fulfilled, 



p E A « 

«e{if'} 



(11) 



7l2 

T12 



cxp 



-0 E A ^ E 



ie{fcf xt } 



7'e{fc! xt } 



exph^fV-Sx)] , (12) 



in which {fcj xt } = 0f xt }u{mf xt } and {fcJf*} = {/| xt }u{m| xt } &ndE x andE Y 
refer to the total internal energy before and after the move, respectively. This 
treatment applies to any simultaneous move of clusters, so that each cluster 
in the decomposition indeed can be moved independently without violating 
detailed balance. It is noteworthy that the probabilities for breaking joint 
bonds in the forward and reverse moves cancel only because the probability 
in the cluster construction factorizes into pairwise probabilities, as opposed to 
the probability for a multiple-particle move in a Metropolis-type algorithm. 

As demonstrated in Fig. 2 for a binary mixture, the results obtained by 
means of the multiple-cluster geometric algorithm agree perfectly with those 
obtained using the single-cluster version. 



4 Performance 

The most striking feature of the generalized geometric cluster algorithm, 
apart from the fact that it creates clusters that can be moved in a rejection- 
free manner, is the speed at which it relaxes size-asymmetric mixtures. We 
illustrate this here for a binary fluid mixture consisting of N s small and N\ 
large spherical particles with size ratio a = cr\/a s > 1. The particles are 
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Fig. 2. Comparison between the single-cluster version (solid lines) and the multiple- 
cluster version (symbols) of the generalized geometric cluster algorithm. The figure 
shows pair correlation functions for the size-asymmetric Lennard- Jones mixture 
described in Ref. [4]. The system contains 800 small (diameter a) and 400 large 
particles (diameter 5a) at a total packing fraction r\ « 0.213. and g^ s represent 
the large-large and large-small correlation functions, respectively. 



contained in a fixed volume, at equal packing fractions j] s = rj\ = 0.1. While 
N\ = 150 is kept fixed, N s increases from 1200 to 506 250 as a is varied 
from 2 to 15. Pairs of small particles and pairs involving a large and a small 
particle act like hard spheres. However, in order to prevent depletion-driven 
aggregation of the large particles, they have a Yukawa repulsion, 



where {3 J = 3.0 and the screening length = a s . In the simulation, the 
exponential tail is cut off at '3<j\. As a measure of efficiency we consider the 
integrated autocorrelation time r obtained from the energy autocorrelation 
function [16], 



For conventional MC calculations, t rapidly increases with increasing a, be- 
cause the large particles tend to get trapped by the small particles. Ac- 
cordingly, an accurate estimate for r could only be obtained for a < 7. By 




(13) 



C(t) = 



(E(0)E(t)) (E(0)) 2 
(E(0y) - (E(0)) 2 



(14) 
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Fig. 3. Efficiency comparison between a conventional local update algorithm (open 
symbols) and the generalized geometric cluster algorithm (closed symbols), for a 
binary mixture (see text) with size ratio a. Whereas the autocorrelation time per 
particle (expressed in /is of CPU time per particle move) rapidly increases with size 
ratio, the GCA features only a weak dependence on a. 



contrast, the generalized GCA has an autocorrelation time that only weakly 
depends on the size ratio, as illustrated in Fig. 3. At a — 7 the resulting 
efficiency gain already amounts to more than three orders of magnitude. 

A crucial limitation of the generalized GCA is that each cluster must only 
occupy a fraction of the entire system. As observed by Dress and Krauth [14], 
the entire system typically will be occupied by a single cluster once the perco- 
lation threshold is reached in the combined system containing a configuration 
and its point-reflected counterpart. For the original system this leads, in three 
dimensions, to a practical upper limit in the packing fraction around 0.23- 
0.25. This number will vary as a function of size asymmetry and, if additional 
pair interactions are present, temperature. It is therefore instructive to study 
the cluster-size distribution as a function of reduced density p* and tem- 
perature T. Figure 4 illustrates the cluster-size distributions as obtained by 
means of the multiple-cluster GCA for a regular (one-component) Lennard- 
Joncs fluid. In the top panel, p* = 0.32 ss p*. Already at temperatures that 
are far above the critical temperature T c « 1.19, the cluster-size distribu- 
tion starts to tend toward a bimodal form, indicative of the formation of 
large clusters. In the vicinity of the critical temperature, the average cluster 
size has become very large. For comparison, the bottom panel displays the 
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Fig. 4. Cluster size distributions as a lunction of relative cluster size X, for a 
monodisperse Lennard- Jones fluid, (a) Reduced density p* = pa 3 = 0.32; (b) re- 
duced density p* — 0.16. Identical symbols compare to identical temperatures in 
both panels. All temperatures are indicated in terms of e/kv- See text for discussion. 



Geometric Cluster Algorithm for Interacting Fluids 



11 



cluster-size distributions at a twice smaller density, p* = 0.16. In this case 
the bimodal shape does not appear until close to T c . 

The properties suggest that the generalized GCA, at least for the one- 
component Lennard- Jones fluid, will not suppress critical slowing - the pri- 
mary advantage of lattice cluster algorithms. For the off-lattice GCA, this 
property is less essential, because of the speed-up it delivers for the simula- 
tion of size-asymmetric fluids over a wide range of temperatures and packing 
fractions. Nevertheless, we have investigated the integrated autocorrelation 
time for the energy at the critical point, as a function of linear system size. 
In Fig. 5 these times are collected for three algorithms. (1) Conventional 
local-update Metropolis algorithm; (2) Wolff version of the GCA; (3a) SW 
version of the GCA, in which each cluster is reflected with a probability 0.50; 
(3b) SW version of the GCA, in which each cluster is reflected with a prob- 
ability 0.75. Just as for spin models, the single-cluster version outperforms 
the Swcndsen-Wang type cluster decompositions. However, all variants of the 
GCA exhibit the same power-law behavior, which outperforms the Metropolis 
algorithm by a factor ~ L 2A . It is important to emphasize that this acceler- 
ation may be due to the suppression of the hydrodynamic slowing down [17] 
caused by the conservation of the density (which may couple to the energy 
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Fig. 5. Energy autocorrelation time r as a function of linear system size for a 
critical Lennard- Jones fluid, in units of particle sweeps, for three different Monte 
Carlo algorithms: Local moves ("Metropolis MC"); GCA with Swendsen-Wang 
type cluster decomposition and probability 0.50 ("SW 50%") and 0.75 ("SW 75%") 
of moving each cluster; single-cluster GCA ( "Wolff" ) . 
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correlations [4]). Another striking point is that already for moderate system 
sizes the generalized GCA outperforms the Metropolis algorithm, despite the 
time-consuming construction of large clusters [cf. Fig. 4(a)] which only lead 
to small configurational changes. 
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Fig. 6. Pair correlation function of dilute colloidal particles (diameter a) in an 
environment of smaller particles (diameter a/5, packing fraction r\ — 0.116) that 
experience a Yukawa-type attractive interaction with the colloids. The symbols 
represent data obtained by means of the hard-core GCA [18]; the solid line was 
obtained from the generalized GCA. 



5 Illustration 



In order to illustrate the capabilities of the generalized GCA, we have com- 
puted the pair correlation function of dilute colloidal particles (packing frac- 
tion rj\ = 0.001) in an environment of particles with a five times smaller 
diameter (packing fraction rj s = 0.116). Large and small particles act as hard 
spheres, but unlike pairs (i.e., large-small) experience a Yukawa attraction 
which promotes the accumulation of small particles around the colloids. This 
system has been studied in Ref. [18] by means of the original hard-core GCA, 
in which clusters are moved according to an acceptance criterion, as proposed 
in Ref. [14]. This potentially greatly deteriorates performance, as entire clus- 
ters will be constructed that are subsequently rejected. Indeed, the authors 
report [18] that the colloid pair correlation function g(r) had to be obtained 
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via numerical differentiation of the integrated pair correlation function, rather 
than through direct sampling. The generalized GCA can handle this system 
without complication, as demonstrated in Fig. 6. In order to obtain data 
comparable to the solid curve in this figure, the authors of Ref. [18] utilized 
a polynomial fit to the integrated pair correlation function, which potentially 
leads to ambiguities. 

6 Conclusion and Outlook 

We have discussed the generalized geometric cluster algorithm, a Monte Carlo 
method for the simulation of fluids by means of geometric operations. This al- 
gorithm creates nonlocal multiple-particle moves that are capable of rapidly 
decorrelating fluid configurations that contain particles of widely different 
sizes. The multiple-particle moves are constructed in such a way that the typ- 
ical decrease in acceptance rate is avoided; every proposed move is accepted 
without violating detailed balance. It is anticipated that this algorithm will 
find widespread application, in particular in the simulation of complex fluids 
and suspensions in which the solvent is modeled as an implicit background. 
Potential generalizations include the treatment of particles with internal de- 
grees of freedom, such as nonspherical particles, and other geometries, such 
as layered systems. 
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